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Abstract 


In this article, we develop an efficient numerical method for one-dimensional 
time-delayed singularly perturbed parabolic problems. The proposed nu- 
merical approach comprises an upwind difference scheme with modified 
graded mesh in the spatial direction and a backward Euler scheme on uni- 
form mesh in the temporal direction. In order to capture the local behavior 
of the solutions, stability and error estimations are obtained with respect to 
the maximum norm. The proposed numerical method converges uniformly 
with first-order up to logarithm in the spatial variable and also first-order in 
the temporal variable. Finally, the outcomes of the numerical experiments 


are included for two test problems to validate the theoretical findings. 
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1 Introduction 


Singularly perturbed problems are widely applicable in different fields of sci- 
ence and engineering. In particular, singularly perturbed convection-diffusion 
parabolic equations, which contain a small parameter and a delay term, are 
effective in the simulation of reservoirs subsurface oil extraction, convective 
heat transfer problems with high péclet numbers, fluid flows, mathematical 
biology, and so on. Problems with delay have also appeared in the study of 
tumor growth, neural networks, and also in problems related to the respira- 
tory system. The dynamics of the delay differential equations may become 
complicated due to time delay as it alters the stability equilibrium of the 
system. In epidemiology, to denote the incubation time or the host infection 
time, delay can be used. From the statistical study of ecological data, it is 
revealed that delay effects are seen in many species’ population dynamics. 
An example of a delay partial differential equation (DPDE) that is used in a 
furnace to process metal sheets is as follows: 
Ow(s, T) 0? W(s,T) 
OT “~~ Gg2 
= v(f(w(s,7 — €))) 


PST) + tg (to(s,7 — §)) ~ (3,7) 

For detailed literature on the above delay model, the reader can refer to [32]. 
DPDEs are closer to real-world phenomena as the solution does not only 
depend on the solution at a current stage but also at an earlier stage. The 
highest order derivative of the singularly perturbed delay differential equa- 
tions (SPPDEs) is multiplied with a small parameter ¢, with at least one 
term included with negative/positive shifts. There are a lot of manuscripts 
available in the literature dealing with DPDEs. For instance, see [14, 29, 22]. 
It is challenging to solve the SPPDEs analytically. Additionally, the pres- 


ence of a small parameter ¢ leads to a sudden change in solution, resulting 
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in boundary and interior layers. Moreover, due to the requirements of an 
unexpectedly large number of mesh points on uniform mesh, a classical nu- 
merical method is not practical to accurately capture the layer in the so- 
lution. In this sense, the scheme mentioned above fails. So, a singularly 
perturbed DPDE numerical solution requires special treatment, which leads 
to the development of a uniform convergent method. Various articles dis- 
cuss the numerical as well as the analytical aspects of singularly perturbed 
differential equations via finite difference and finite element methods. For 
example, see [2, 7, 19, 28, 24, 4, 23, 25]. However, only a few articles are 
available on partial differential equations with delay arguments in which an- 
alytical and numerical solutions are discussed. Ansari, Bakr, and Shishkin 
[1] used a finite difference method (FDM). Their method is implicit in time 
and centered in space on a piecewise uniform mesh. They achieved an order 
of convergence O(( N~'InN)? + M~*) for the singularly perturbed delay 
parabolic reaction-diffusion problem. Gowrisankar and Natesan [10] solved 
SPPDEs of a convection-diffusion type by virtue of the FDM. In this ar- 
ticle, the spatial domain through the adaptive mesh and temporal variable 
using implicit Euler has been discretized, and the method converges with or- 
der O((N~'In N) + At). Gowrisankar and Natesan [21] established uniform 
convergence for a singularly perturbed convection-diffusion problem with a 
convergence of order O(At + N~!*7), where 0 < g < 1. Kumar and Ku- 
mari [15] constructed an ¢-uniform convergent FDM for the reaction-diffusion 
initial-boundary value problem with a delay in the temporal direction. The 
authors have confirmed the order of convergence for a fully discretized scheme 
is O(At + N-? In® N) utilizing extended cubic B-spline on the layer-adapted 
mesh. Kumar et al. [16] have explored uniform convergence of time-delayed 
parabolic partial differential equations (PDEs) on graded mesh generation 
algorithms based on entropy function. Das, Govindarao, and Mohapatra 
[3] proposed a weighted monotone numerical scheme comprised of Crank- 
Nicolson in time and a weighted monotone hybrid scheme for spatial deriva- 
tives on Shishkin mesh, which is parameter uniform converges with order 
O((N~tIn N)? + (At)?). Gartland [8] studied boundary value problems on 
the graded mesh in which the scheme converges of order O(h*), where k 


times as many points are available inside the layer as outside. The con- 
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vergence of order two in both the spatial and temporal variables has been 
established by Gupta, Kadalbajoo, and Dubey [11]. The difference scheme 
on the Bakhvalov-Shishkin mesh is used in [12] to explore the first-order con- 
vergence. Indeed in all these works, authors have described delay /nondelay 
problems on layer-adapted mesh only. There are various numerical meth- 
ods used in different branches of science and engineering; for reference, see 
(26, 27, 5, 30, 31, 6, 33). Any work related to the convergence of difference 
schemes on the modified graded mesh has not been noticed yet. Therefore, we 
are now in a position to establish a different scheme on the modified graded 
mesh. Motivated by the work of [1, 10, 4, 23], we proposed an upwind FDM 
on a modified graded mesh for a time-delayed convection-diffusion parabolic 
problem. We also establish that the method is almost first-order convergence 


without a logarithmic factor in the mesh parameter NV. 


Many articles discuss convection-diffusion problems involving bound- 
ary layers with layer adaptive meshes like Shishkin, Bekhvalov, and adap- 
tive meshes. Gowrisankar and Natesan [10] solved SPPDEs of convection- 
diffusion type by the FDM on piecewise uniform mesh, which converges 
uniformly with an order O((N~!In N) + At). Kumar et al. [17] analyzed 
convection-diffusion problems on the adaptive mesh via equidistribution of 
monitor function and obtained first-order convergence. Generating adaptive 
meshes for the parabolic problems requires the iterative process to obtain the 
grid at each time level. So, it requires high computation in comparison to 
other meshes discussed above. It is clear from the numerical findings that 
the modified graded mesh proposed in this work provides the equivalent rate 
of convergence as the adaptive meshes with negligible computation in com- 
parison with the adaptive algorithm. Also, the modified graded mesh helps 
improve the accuracy and efficiency of the method. 

The description of the problem is as follows: Let 8 = A x (0,7], where 
A = (0,1) and Y = TY; U YT, U T,., in which YT, and TY; are, respectively, the 
right and left side of the rectangular domain X and YT; = (0, 1] x [—€,0]. The 
following problem on 8 = A x (0,7) is examined in this paper: 


—d(s,T)W(s,7—€)+g(s,7), forall (s,7)EX, (1) 
(s,r) € {0} x TY; = {(0,7),0<7< T}, 


tS 
a 
= 
ers, 
& 
4 
ls 
lI 
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w(s,T)=W,(7T), wW(s,7)€ {1} x YT, ={(1,7),0<7 <TH, 
w(s,7) = WV,(s,7), (5,7) € Ty = [0,1] x [-€, 0], 


where 


LtW(s,7) = (2 rele ee b(s,7)) ws, r). (2) 

Or Os? Os 
The above-considered singular perturbation parameter is ¢ € (0, 1]. Addition- 
ally, € > 0, a delay parameter, satisfies the equation T = gé for some positive 
integer q. Here, the above considered functions &(s,7), 6(s,7), d(s,7), 9(s,7), 
W,, U,, Uy, and Wy are sufficiently smooth, bounded and also fulfill the fol- 


lowing conditions: 
a(s,T) >6>0, b(s,T) >0, and d(s,T) >B>0. 


The reduced problem corresponding to (1) is 


eee + a(s,7) HE + b(s,7)t0(s,7) 

= —d(x,7)tio(s,7 — €) + 9(s,7), for all (s,7) EX, 
wo(s,7) = Va(s,7), (s,7) € Ty = [0,1] x [-€, 0], 
Wo(s,T) = U,-(s), (s,T) € T>. 


(3) 
By considering s to be constant, the solution to the reduced problem (3) is 
vertical lines, which ensures that the boundary layer that arises will be of the 
parabolic type. The compatibility condition for the initial functions V,(s,7) 
are also satisfied at the corner points (0,0), (1,0), (0, —€), and (1, —), 


(0,0) = W(0), 
W, (1,0) 7, W,(0), 


d¥;(0) _82W,(0,0) 
dt . Os? 


= —d(0,0)¥,(0, —€) + g(0,0), 
cso pe, ow, (1, 0) + b(1, 0) Wo (1, 0) 


= —d(1,0)%,(1, -€) + g(1,0). 


+ 5(0, 0) Wy (0, 0) 
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Notation: Assume that C is used as a generic constant throughout the 
paper. Also, C is independent of the perturbation parameter ¢ and the mesh 
point N. We use discrete maximum norm to study convergence, which is 


defined as 


I|\|a = max |t0(s)|- 


The Paper is presented as follows. Continuous solutions and their derivatives 
are covered in Section 2. The constructed numerical method is presented 
in Section 3. Section 4 studies the error analysis of the proposed method. 
Numerical examples and results support the theoretical findings in Section 


5. The paper ends with a conclusion presented in Section 6. 


2 Bounds for continuous solution and its derivatives 


This section discusses the bounds for the solution and its partial derivatives 
across the prescribed domain. Furthermore, for the proof of ¢-uniform con- 
vergence in a later section, we split the analytical solution w of (1). Also, 
it establishes strong bounds on the layer and smooth component, which we 
obtain by splitting the analytical solution. The operator L7 defined in (2) 


satisfies the following maximum principle. 

Lemma 1 (Maximum principle). [15] Assume that €(s,7) > 0 holds true for 
all values of (s,7) € Y. If L7&(s,7) > 0 for all (s,7) € X, then €(s,7) > 0 for 
all (s,T) ER. 


Proof. Suppose (s*,7*) € ® such that 
€(s",7") = min &(s,7), 


and assume £(s*,7*) <0. Also, we have d6(s" sr") = 0, O6(8" sr") = 0, and 
oes") > 0. Then 


LyE(s", 7") <0, 


which is contrary to our supposition. It follows that €(s*,r*) > 0 and so 
€(s,T7) > 0 for all (s,7) ER. 


Lemma 2 (Uniform stability bounds). [13, 1] For the specified function 7 in 


the defined domain of the differential operator LZ in (2), we have 
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IIn|| < A+ eT) max{||Lsnll, lll, ll, loll 
and for any solution w(s,7) of (1), we have 
I|(s,7)|| < A + oT) maxfllgll, Mall, rll, Wold, 
where a = max{0,1— 6}. 


Proof. We define barrier function as 


o(s,7) = (1+ aT) max{|lgl}, [Y(7)I, rr) Mor) + Ws, 7), 
for all (s,T) ER. 


Then, at initial condition 


nt 
= 


wy (s,0) = (1+ aT’) max{|lgl|, || ¥1(0) |], | Y-(O) |], WoO) |} + w(s, 0), 


we obtain 


w~(s,0) >0, forall (s,r) € To. 


At boundary points, we have 


y~ (0,7) = 1 + aT) max{|[gl], Yer) I], Mr) IL, Mo) + WO, 7), 


which gives 


w~ (0,7) >0, for all (s,7) € TY). 


Similarly, we have 


w~(1,7) > 0, for all (s,7) € T;. 


The differential operator L7 satisfies 


LUw* (8,7) = (Wr(8,7) — e(+Ws5(s,7)) + A(s,7)(4ts(s,7)), 

+ b(s,7)((1 + aT) max {ll [| ¥i(7)II, Yr) I Wo(7) I} + BCs, 7))) 
= (40, (s,7)Fe(tss(s,T))+4(s, T)t5(s, 7) 46(s, T)W(s, 7), 
+ b(s,7)((1 + aT) max {|| gl, [| ¥i(7)II, Yr), I Wo(7) I} + BCs, 7))) 
= LZ W(s,7) + 6(s,7)((1 + aT) max |lgll, |Ye(7)I, Ir IL, IoD 
> 0. 


Thus, the maximum principle asserts that 


Iran. J. Numer. Anal. Optim., Vol. 14, No. 1, 2024, pp 77-106 


Kumar and Gowrisankar 84 


w*(s,T)>0, forall (s,r) ER. 
Finally, we get 


I|(s,7)|| < A + aT) maxfl|gll, Mall, Yell, Wold. 


Theorem 1. Let the data a € C@+11+7/2)(A), 6, d, g € C2+71147/2) (R), 

Ge OP e(0,7), we C707) eee ey). S (0.1), 

meet the necessary compatibility requirement on the corner. Consequently, 

(1) has a unique solution # and w € C(4+1?+7/2) (&). Additionally, bounds 

on derivatives of solution w(s,7) is satisfied for all nonnegative integers i and 
j, where 0 <i 427 <4, 7 

Oiti yw 

laraal| sce 

Proof. The proof of existence and uniqueness can be found in [18]. To estab- 


lish the derivative bound, the variable s is replaced to ev = s and follows the 


approach given in [1]. 


In the above Theorem 1, we establish the bound for exact solution w(s,7), 
which is not sufficient to establish the c-uniform convergence of the method. 
In order to get strong bounds on the derivatives of the exact solution w, we 


decompose the above w into smooth and layer component, 


w(s,T) = 0(s,7) +w(s,T), (8,7) ER, 


where the following differential equations are satisfied by the smooth compo- 
nent, 6(s,7), 


L76(s,7) = —d(s,7)d(s,7 —€) + g(s,7), (5,7) EX, 


with initial condition 


and the boundary condition 


6(0,7) = w(0,7), 61,7) = W017), O<7T<T. 
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The smooth component 6(s,7) is further decomposed as 


0(s, 7) = b0(s, T) + 61 (s,T), (s,T) ER, 


where 69(s,7) is the solution of reduced problem 


ols-7) | B(s,7)00(8,7) 


= —d(s,7)6o(s,T ~~ f) a G(s, #), for all (s, T) € X, 


= —d(s,r)61(s,7 é)+ O°éo(s,7) 


o1(s,7T) =0, (8,7) € To. 
Furthermore, 0 satisfies 


L76(s,7) = —d(s,7)0(s,7 —7T)+4(s,7), (s,T) ER, 
6(s,T) = W(s,7), (s,T) € To, 
6(0,7) = 60(0,7), 6(1,7) = 6o(1,7) + e61 (1,7), T'S (0, T], 


and the singular components w satisfies 


L’w(s,7) = —dw(s,r — €), (8,7) ER, 
; (s,T) € To, 
w(0,7) = 0, w(1,7) = w(1,7) — o(1,7), 7 € [0,7], 


and the estimation of singular and the smooth component are discussed in 


Theorem 2 given below. 


Theorem 2. Assume that & € C(4+7-2+7/2) (A), 6, d,g € C4+12+7/2) (R), Wy € 
C8+9/2((0,T]), v, € C7 ((0, T]), UV, € Clo+13+7/2)(T,), y € (0,1), fulfill 
the necessary compatibility condition on the corner. Then, for all nonnega- 
tive integers i,7 such that 0 < i+ 27 < 4, we have the following bounds for 


the smooth 0 and the layer part w in the decomposition of solution w is 
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OrdG : 
1 1-7 
laze < Sou re") 
and ae 
a Jw 
< -i —ys/E 
locas Os'Ord | < Cee 


Proof. For the proof of theorem, one may refer [20]. 


3 Grid construction and numerical discretization 


In this section, grids are developed to discretize the problem in both spatial 
and temporal directions. Later, we use a backward Euler in the temporal di- 
rection and an upwind method on the spatial derivative for the discretization 
of the problem (1). 


3.1 Temporal discretization 


To establish the convergence of (1) at each instance, we use the uniform time 
grid. 
OM ={,=kAr, k=0,1, ..., M, ArM =T}. 


Here M represents the grid points in the temporal direction. 


3.2 Spatial discretization 


We generate a modified graded mesh, AN in the interval [0,1] as follows: 


Ho = 0, 
by = 24, 1 

Mju =uj(1+ph), F<j<Nn-2, 
un =1, 


where the parameter h satisfies the following nonlinear equation: 
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In(1/e) = (N'/2) In(1 + ph). (5) 


The above selection of the parameter h ensures that there are N/2 grid 
points in the interval [e,1], which are distributed gradely in the interval 
[e, 1]. Numerical verification stimulates that in comparison to ({4n—2, UN-1), 
the interval (~y—1,1) is not too small. In the subinterval (0, <], we distribute 
N/2 points with uniform step length 2¢/N, while in the subinterval [e, 1], 
we first find h for some fix N by means of the nonlinear equation (5), and 
corresponding to that h, we distribute N/2 points in [e, 1]. 

The mesh length is denoted by h; = uj; — wj—1 for 7 = 1,2,...,N. 


Jah ++} +j—}-_+- + + ____| 


tt) é 1 


joa te -N/2 ‘| 


Figure 1: Modified graded mesh for N = 32 and ¢ = 107}. 
Remark 1. The mesh size in piecewise uniform and the modified graded 
region is given by 


2e/N for j=1, 2, ..., N/2, 
phpj_1 forj=N/2+1, N/2+2, ..., N. 


Lemma 3. The mesh defined in (4) satisfies the following estimates: 


0 for j =1, 2, ..., N/2, 
|Ajta — byl S 
Ch for j=N/24+1, N/2+2, ..., N. 
Proof. Initially, we consider 7 = 1, 2, ..., N/2. As the mesh is uniform in 


this portion, so nothing to prove. 
For j7 = N/2+1, N/2+2, ..., N. We have 


|nj41 — hy| = |phy; — phyj—s| 


= ph| pj — Mj-1| 
= ph? wj-1 
< Ch. 


Here, we have taken 0 < p,h < 1. 
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Lemma 4. For the modified graded mesh defined in (4), the parameter h 


satisfies the following bound, 
h< CN In(1/e). 


Proof. Let K, denote the number of partitioned points in (4) such that yp; < 
é, for j = 1, 2, ..., N/2. Clearly, Ki < C/h. Let Ko be the number of 
points in the partition (4) such that yj; > ¢. Let yy241 be the smallest 
point such that y; > ¢. We have to estimate the bound for Ky. Assuming 
ph <1, we have 


N/24+1 N/24+1 j 
N 1 Mj+1 

= yt fa 

N/24+1 Hj 


lI 
> 
aS 
= 
a. 
| 
ws 
ae 
Q. 
= 


IA 
bo 
™~ 
> 
= 
= 
a. 
S 
aa 
a 
i 
= 


< 2(ph)-* In(1/e). 
Recalling N = Kk, + Ke, we have 


N <C/ph+2(ph)~" In(1/e) 
N < 1/h(pC + 2pIn(1/e)) 
N < 1/h(C In(1/e)). 


Finally, we get 
h<CN'In(1/e), 
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where N represents grid points in the s-direction. 


Now, we systematically discretize our domain. For this, a non-uniform 
grid AX on A with N mesh points is obtained by distributing N/2 points 
uniformly in the layer part and N/2 points nonuniformly outside the layer 
part. Here, we consider an equidistant grid Q” and ©! with uniform step- 
length Ar on [0,7] and [—€, 0] with M and / grid points, respectively. Then, 


the discretized domain will define as 
RV = AN x OM, TR=AN xl, 


and the boundary points TX of 8X are YY = 8X 9 Y. Additionally, es 
= Y, and TX = x A YT, define the left and right boundary points, 
respectively. Also, NN = AN x OQ, where Q!, is the set of uniform mesh 
in [(i — 1)&, 7]. Now, we discretize our problem with the above-discretized 
mesh. We apply an upwind difference method and backward Euler for spatial 
and temporal derivatives. Then 


J 


a ee a . 
(Dz - 67272 + (j,k +1) DES; 41 
+6;,64155,h41 = (j,k + 1)Sjn—-1 + 95,6415 
So,et1 = Vi(te+1), (6) 


Swirti = Ur(Tr41), 


Sip = Va(8;,—Tp), j=l, 2, ..., N-1, p=0,1, ..., 1. 


By arranging (6), we will get a tri-diagonal system 


(05,441) S5—1,k+1 + (Bj,6+1)S5,h41 + (5,641) Si4tkt1 = j,k, (7) 
where 
—eA 
Aj ko = —, 
Tk+1 hgh, 
eAtT eAr a(j,k +1)Ar 
byez =14 ae ea U ) b5,n41 AT, 
hy4ih; hh; hj+1 
—eAr a(j,k +1)Ar 
Vik+1 = al ) 
hjzih; Aj+a 


hyn = —d(j,k +1)S5 4 -1AT + 95,41 AT + Six, 
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2 ie 
hj = + a 


and the step length Az satisfies JAr = €, and also Tt, = pAt, p > —l, and | 


is a positive integer and the mesh function U(s;,7%) = U;,x is defined by 


Uj+1,k — Uj,k Uj,k — Uj-1,k 


a ee. a 
De U3,b = ys » DS Uj.KR = his ? 
g@+1 | 
= Uj, — Uj,b—-1 
eas ied Js 
D; Uj,k —_ h; ) 
g+1 
and 2 
520. _ (D3 ae D; )e; k 
sUj,k = x 
h 


3.3 Numerical algorithm 
The following algorithm provides the grid construction and the corresponding 
numerical solution: 


Step 1. Given the number of mesh points in the temporal and the spatial direc- 
tion, M and N, respectively, take uniform mesh points in the temporal 


direction, {7% }}“o- 


Step 2. For the finer part in the spatial direction (i.e., [0,¢]), we have the uni- 
form mesh te 


Step 3. For the coarser part (i.e., [¢,1]), the graded mesh parameter h is ob- 


tained by solving the nonlinear equation (5) by the bisection method. 


Step 4. Using the graded mesh parameter h, obtain the graded mesh in the 


interval [e, 1] from (4). 
Step 5. Set k = 1. 


Step 6. For the value of k, solve the tridiagonal system (7) to obtain the solution 
for the time level t = k. 


Step 7. k =k-+1 goto Step 6. 


Step 8. If k = M, then stop and mark S;,%, as the required solution. 
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4 Error analysis 


Lemma 5 (Discrete maximum principle). [9] Assume that the mesh function 
(8;, 7k) satisfies €(s;,7,) > 0 on YY. If LIN E(s;,7%) > 0 on (83,7) ERY, 
then €(s;,T) > 0 on an 


Lemma 6. For any solution S(s;,7,) of (6), we have 
||3(sj, 7) |] < (1+ aT) max(||L5™ |], || Yl). (8) 
Proof. By constructing the barrier function, we have 


S*(sj,Th) = (1+ @2) max(||LP" ||], ||Wllew) + $(s;,7%)- (9) 


We can obtain (8) by applying (5) on (9). 


Theorem 3. Let w# and S be the solution of continuous problem (1) and 
discretized problem (6), respectively. Furthermore, the both solutions meet 


the corners compatibility requirements. Then, the error estimate is given by 
max |(t — S)(s;,T)| < C[AT + N7*In(1/e)], (sj, Th) ERY. (10) 
Here, constant C is free from N, Ar, and e. 


Proof. The proof is similar to the proof given in [10, 1]. We derive it by using 
our mesh briefed in sections 3.1 and 3.2 and (5). 
To prove the theorem for different time levels, we first divide the domain & 
into X = N; UNe, where X; = A x [0,€] and Ny = A x [€, 2€]. The discretized 
domain R™ = RY «Mo. where ST = AT xO and RF = APY xo, 
First, for r € (0, €], the right-hand side of (1) becomes —d(s,r)(s,r — €) + 
g(s,7), which is free from ¢. Thus the results in [10] is applicable, and we 
obtain 

max |(w# — S)(s;,7%)| < C[Ar + N~' In(1/e)]. (11) 


For t > &, the term w(s,7 — €) is not free from ¢. So, we are required to 
examine the estimates between the numerical solution S and the analytical 


solution w over the interval [€, 2]. The following SPPDEs is considered: 


a) oo. Os» . 
(5 en. G(s,7) 5 + b(s,7)) W(s,7) 
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= —d(s,T)W(s,7 —€)+ (8,7), (8,7) € ny, (12) 
w(s,7) = w(s,71), 5 € (0,1), 


w(0,t) = Wo(r), wWil,r7)=Vi(r), 7 € [€, 26]. 


To determine the numerical solution, we discretize (12) by using upwind finite 


difference for spatial derivatives and the backward-Euler for time derivative, 


L7S(83,Tk) = D7 Sj% — £62554 + 4; 4D, S35, k + by 055, 
=— 1 42555, —1 Ol Site “C85 7K) Se cre (13) 
S(8;,Tk) = $1(8;,Tk), (8;,Tk) € n, 
S(0, 7%) = Wol(tTr), 
S(1,7%) = Wilre), Te € 4, 


where $;(s;,7,) is the approximate solution obtained over the interval Ny’. 
Now, we divide the solution w of (1) into its regular and singular components 
as W = 4+ 2, and furthermore, 9 = yo + €y1, where yo is the solution to the 


reduced problem, 


(2 + a(s,7)2 + i(s,7))yo(s,7) 
= —dyo(s,7 —€)+9(8,7), (8,7) € (0,1) x (6,28), 
yo(s,7) = W(s,7), (8,7) € (0,1) x [0,4], 
yo(0,7) = 00,7), 7 € [€, 26], 


and 


. oO? : 
Lryi(s,7) = —d(s,T)yi(s,7 —€) + ils =) 


yi(s,7) = 0, (s,7) € (0, 1) x (0, €], 
yi(0,7) =yi(1,7) = 9, TE [E, 2€]. 


Furthermore, % satisfies 


L59(3,7) _ —d(s,7)9(s,7 = f) + Olay), (s,T) € (0, 1) x fest 
9(s,T) = w(s,7), (8,7) € (0, 1) x [0, €], 
y(0,7) = yo(0,T), g(1,7) = yo(1,7), TE [f, 2€]. 
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The singular components z satisfies 
£5 2(8,7 £) = 2(8,7 f)3 (8,7) € (0, 1) x (ear 
2(s,7)=0, (8,7) € (0,1) x [0,€], 
0, 2(0,7) = Wi(r)—yo(0,7), 7 € [€, 26]. 


Now, we decompose the numerical solution S of (13) into S = Y+Z. Here, Z 
represents the singular component of the decomposition, and Y, the regular 


component, is the solution to the following non-homogeneous problem: 


hae = —dY (s;, Tk-1) + g(34, Tk)s (53,7) € ny, 
Y (85, Tx) = 51 (83, Tk), (8;, Tk) € (0, 1) x (0, ), 
Y (0, Tr) = g(0, Tr), Y (1, Tr) = 9(1, Tr), Tk © eA 


and the singular component Z must satisfy 


ied = —dZ (83,71), (8;, Tr) € ny, 
VACIRED = 0, (8457) S ny, 
Z(1, Tr) = 0, Z(0, Tr) = Wi (Tk) _ g(0, Tr), Tk € 08 


Therefore, the error can be written in the form 
S-w=(Y-9)+(Z-2). 


Now, we will establish the bounds for the smooth and the layer components. 
To establish the bound for the smooth component, we use the classical argu- 


ment. The smooth error component can be expressed as 


LUN (Y — 9) = —dY (83, Te-1) + 9(83, 76) — LENG 


= —dY (s;, Tk—1) +r Ly + dg(s;, Tk—1) _ LING 
= d( (83, Te—1) — ¥ (85, 7-1) + LONG. 


Thus we obtain 


te Os 

LEN (Y — 9) = d(G(sj, 71) — ¥ (83, 74-1)) — (aq — 85)9 
Oe he ae 
1 (5. _ 67)9 an (83, Tk) (a— _ 05 )g. 
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When absolute values are taken on both sides and used with (10), the above 


inequality becomes 


; - o ; 
ILEN(¥ — §)| < C(Ar +. N7Mn(1/e)) + el( 55 - 82)8 
) 


= 0 ; 
+5 -85)0 + (== 65)0 


Using the Taylor series expansion, we have 


T yy : 2 
JES" (¥ — 9) < C(Ar +N In(1/e)) + (8541 — 1-1) | 55a 


?G)|_. (te —Te-1) || PG 
32 || : 


= vol 2 lldr 


On using the estimates of derivatives, bounds of mesh length, and discrete 


maximum principle, we have 
\(Y — 9)(8;,7%)| < C[Ar + N~* In(1/e)]. (14) 


Like continuous z, Z is discretized for estimating the singular component. 
Thus 


Lg => —dZ(s;,Tr-1); (35, 7c) E ne, 
Z(8;,Tk) =0, (85,TR) € Re 


Z(1, Tr) =0, Z(0, Tr) _ Wi (Tk) = 9 (0,7), Tk E Ob: 
The singular component error can be expressed as 


LON 22) LO" 2 LONe 
= —dZ(si,Tn—1) - LE" 3 
= (LE -L3")2 


Then, the classical estimates gives 
heal ea £)(8;,Tk)| < C[Ar + N7* In(1/e)]. 


The discrete maximum principle is satisfied by the operator L7 and also, 


due to the uniform boundedness of the inverse operator, the above inequality 
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reduces to 
\(Z — 2)(83,T)| < C[Ar + N71 In(1/e)]. (15) 


Combining (14) and (15) completes the proof on [€, 2€]. In the similar fashion, 


we estimate the error in the successive interval in time by using induction. 


5 The numerical result with examples 


This section consists of two test examples with boundary layers to illustrate 
the numerical method discussed above. Using tables and graphs, we present 
the findings of numerical methods. All the results are performed by taking 


p =0.9. Numerical results confirm our theoretical findings. 


Example 1. Consider the following problem: 


= W0(s,7—1)+9(s,7), (s,7) € (0,1) x (0, 2], 
w(s,0) = wo(s,T), (s,7) € (0,1) x [-1, 0], 
w(0,7) =7T, w1,r)=0, 7 € (0,2). 


etn +2rs cos((7s)/2), we will 


Using the exact solution, #(s,7) = G—==t.— 


obtain wo(s,7) and g(s,7). Additionally, the maximum point-wise error for 


each € is defined by 


eg’? = max |(t — $)(sj,7K)| (83,7) ERY, 
where w and S are exact and approximate solution, respectively. The order 


of convergence is computed by 


NAr loge’ St feet) 
P = 


log 2 


For Example 1, the computed maximum pointwise error in the layer region 
is shown in Table 1, and the associated convergence order is in Table 2. 
Additionally, in the smooth region, Table 3 shows the maximum pointwise 
error, and Table 4 presents the related order of convergence. The plots for 


the numerical solution of Example 1 for the parameters N = 128, « = 107? 
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and ¢ = 107° are shown in Figures 2a and 2b. At s = 0, the figure shows 
the existence of a boundary layer. Error plot of Example 1 for N = 128, 
e = 10-2, and e = 10~® is presented in Figure 3a and 3b, respectively. 
Figure 6a is a plot of Example 1 on the log-log scale for maximum pointwise 
error. From Table 5 and Figure 6a, it is clear that the applied scheme is 


first-order uniformly convergent. 


Table 1: Maximum error on the modified graded mesh in the layer region for Example 
1 


Number of intervals N/Ar 


é 128/45 256/45 512/45 1024/3, 2048/si5 4096/ 355 


10-? 2.2150E-01 1.0252E-01 4.8929E-02 2.3896E-02 1.1828E-02 5.8921E-03 
10-4 5.0063E-01 2.3021E-01 1.0696E-01 5.1049E-02 2.4870E-02 1.2266E-02 
10-6 7.6497E-01 3.6319E-01 1.6700E-01 7.8479E-02 3.7821E-02 1.8535E-02 


10-8 9.9761E-01 4.9961E-01 2.3016E-01 1.0699E-01 5.1072E-02 2.4881E-02 


Table 2: Convergence rate for Example 1 using a modified graded mesh in the layer 


region 
Number of Intervals N/Ar 
€ 128/35 256/35 512/45 1024/5 2048/75 
10° 1.1158 1.0731 1.0373 1.0158 1.0053 
10-4 1.1057 1.1067 1.0703 1.0394 1.0207 
10-6 1.0447 1.1156 1.0919 1.0556 1.0302 
10-8 0.9509 1.1049 1.1060 1.0697 1.0391 
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Table 3: Maximum error on the modified graded mesh in the smooth region for Example 
cy 


Number of Intervals N/Ar 


é 128/45 256/35 512/45 1024/4 2048/ za5 4096 / 345 


107? 2.2150E-01 1.0252E-01 4.8929E-02 2.3896E-02 1.1828E-02 5.8921E-03 
10-4 5.0063E-01 2.3021E-01 1.0696E-01 5.1049E-02 2.4870E-02 1.2266E-02 
10-° 7.6497E-01 3.6319E-01 1.6700E-01 7.8479E-02 3.7821E-02 1.8535E-02 


10-® 9.9761E-01 4.9961E-01 2.3016E-01 1.0699E-01 5.1072E-02 2.4881E-02 


Table 4: Convergence rate for Example 1 using a modified graded mesh in the smooth 


region 
Number of Intervals N/Ar 
€ 128/45 256/35 512/45 1024/2, 2048 / 755 
10-? 1.1113 1.0672 1.0339 1.0146 1.0053 
10-4 1.1208 1.1058 1.0671 1.0374 1.0198 
10-° 1.0747 1.1209 1.0894 1.0531 1.0289 
10-8 0.9977 1.1182 1.1051 1.0669 1.0375 


(a) N = 128, e= 10~? (b) N = 128, e= 10-8 


Figure 2: Solution profile for Example 1 for various values of N and ¢ 
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Table 5: Maximum error and rate of convergence on a modified graded mesh for Example 


1 


Number of Intervals N/Ar 


€ 128/75 256/35 512/75 1024/5, 2048/3245 4096/3355 

0-2 2.2150E-01 —-1.0252E-01  4.8929E-02--2.3896E-02 :1.1828E-02_5.8921E-03 
1113 1.0672 1.0339 0146 1.0053 

0-4 -5,0063E-01 -2.3021E-01 —-1.0696E-01_-5.1049E-02 -2.4870H-02 —-1.2266E-02 
1208 1.1058 1.0671 0374 1.0198 

0-® —7.6497E-01 —3.6319E-01 6700E-01 7.8479E-02 3.7821E-02 —-1.8535E-02 
0747 1.1209 1.0894 0531 1.0289 

0-8 9.9761E-01  4.9961E-01 -2.3016E-01_—«1.0699E-01 5.1072E-02-—-2.4881E-02 
0.9977 1.1182 1.1051 0669 1.0375 
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Error 


0.5 


u t 


(a) N = 128, e= 107? (b) N = 128, e= 1078 


Figure 3: Error profile for Example 1 for various values of N and ¢ 


Example 2. Consider the following problem: 


(8,7) — eWss(s,7) — 200.(s, 7) 

= t(s,7 — 1) +107? exp(—7)s(1— 5), (8,7) € (0,1) x (0,2), 
w(s,7) =0, (s,7) € (0,1) x [-1, 0], 
w(0,7) =0, w(1,7)=0, 7 € [0,2]. 


In Example 2, we do not have an exact solution. To determine the max- 
imum pointwise error and to establish an order of convergence, we employ 
the double mesh principle to the computed solution. For this, we consider 
2M and 2N mesh interval in the temporal and spatial directions to deter- 
mine the numerical solution $(s;,7,) on the mesh X2N = A2% x 02™, for 
j=1, 2, ..., N. Additionally, the jth mesh point of AN and the 2jth mesh 
point of A2 are identical. Furthermore, for each ¢, the maximum pointwise 


error is defined by 
ENAT = max |(S— 5)(sj,7)|, (8j.7%) ER, 
and the order of convergence is given by 


sini: loge"? fer er7?) 


log 2 
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In Example 2, the computed maximum pointwise error is shown in Table 
6, and the associated convergence order is in Table 7. The plots for the 
numerical solution of Example 2 for the parameters N = 128, « = 107? and 
e = 10-8 are shown in Figure 4a and Figure 4b. The boundary layer at 
s = 0 can be seen in these figures. Error plot of Example 2 for N = 128, 
¢ = 10-7 and ¢ = 10-® is presented in Figure 5a and Figure 5b respectively. 
Figure 6b is a plot of Example 2 on the log-log scale for maximum pointwise 
error. Table 7 and Figure 6b show that the employed scheme is first-order 


uniformly convergent. 


Table 6: Maximum errors for Example 2 on modified graded mesh 


Number of Intervals N/Ar 


é 128/45 256/35 512/% 1024/4 2048/45 4096/55 


10-2 1.4063E-03  6.3978E-04  2.9302E-04 1.3810E-04 6.671 7E-05 3.2746E-05 
10-4 2.8608E-03 1.4120E-03 6.4327E-04 = 2.9485E-04 1.3905E-04 — 6.7192E-05 
10~® 3.8686E-03 2.1730E-03 1.0218E-03 4.6404E-04 2.1513E-04 1.0258E-04 


10-8 4.3736E-03 2.8588E-03 1.4119E-03 6.4336E-04 2.9491E-04 1.3906E-04 


Table 7: Rate of convergence for Example 2 on a modified graded mesh 


Number of Intervals N/Ar 


€ 128/4 256/35 512/45 1024/3, 2048/45 
10-? 1.1362 1.1266 1.0853 1.0496 1.0267 
10-4 1.0187 1.1343 1.1254 1.0844 1.0493 
10-6 0.8321 1.0885 1.1389 1.1090 1.0685 
10-8 0.6134 1.0178 1.1340 1.1253 1.0846 


6 Conclusions 


This article explored the upwind difference scheme on a modified graded 


mesh. Optimal error bounds were established in the maximum norm. The 
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(a) N = 128, e= 107? (b) N = 128, e= 1078 


Figure 4: Solution profile for Example 2 for various values of N and ¢ 


x10 


Error 


0.5 


u t 


(a) N = 128, e= 107? (b) N = 128, e= 1078 


Figure 5: Error profile for Example 2 for various values of N and ¢ 


method was proven to converge with first order up to logarithm. It can 
be observed from the numerical simulation that this logarithm term In(1/e) 
is not significant. Analytical estimations were carried out to gain uniform 
convergence results of computed solutions in the article. Two test problems 
were included to check the efficiency of the numerical method. Tables have 
been used to present the maximum pointwise error and order of convergence. 


Finally, numerical experiments confirmed the theoretical findings. 
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—o-e=107? 
—+—e=104 
—A~-e=10% 
7 —<—<=10% 
SN | 0 n(1/2)) 


Max-Error 
a" 
S 
Max-Error 


10-8 10° 
2 3 4 
10° 10 10 102 10° 


N N 


a) Log-log plot for Example 1 b) Log-log plot for Example 2 


Figure 6: Log-log plot for Examples 1 and 2 
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